Analysis of fractional order model on higher institution students’ anxiety towards mathematics with optimal control theory

Anxiety towards mathematics is the most common problem throughout nations in the world. In this study, we have mainly formulated and analyzed a Caputo fractional order mathematical model with optimal control strategies on higher institution students’ anxiety towards mathematics. The non-negativity and boundedness of the fractional order dynamical system solutions have been analysed. Both the anxiety-free and anxiety endemic equilibrium points of the Caputo fractional order model are found, and the local stability analysis of the anxiety-free and anxiety endemic equilibrium points are examined. Conditions for Caputo fractional order model backward bifurcation are analyzed whenever the anxiety effective reproduction number is less than one. We have shown the global asymptotic stability of the endemic equilibrium point. Moreover, we have carried out the optimal control strategy analysis of the fractional order model. Eventually, we have established the analytical results through numerical simulations to investigate the memory effect of the fractional order derivative approach, the behavior of the model solutions and the effects of parameters on the students anxiety towards mathematics in the community. Protection and treatment of anxiety infectious students have fundamental roles to minimize and possibly to eradicate mathematics anxiety from the higher institutions.

develop negative attitudes and emotions toward mathematics that will consequently lead to a decreased level of achievement in mathematics 4 .
Studies of real world situations in different disciplines of social sciences, natural sciences, engineering and technology, businesses, and even the arts using mathematical modelling approaches have attracted the attention of various scholars because they provide a better understanding of their investigation, existence, stabilities, and controlling mechanisms 3,10,11 . In the past, a lot of mathematical models of real world situations are consist of a system of integer-order differential equations. Nevertheless, in the last few decades, the fractional-order differential equation approach has been applied in the modelling of real world situations since this approach provides a greater degree of accuracy than the integer order differential equation approach 10,11 . Recently, researchers' interest to conduct a fractional order modelling research to investigate real world situations phenomenon in different disciplines increases because of a lot of properties which are not found in classical integer order modelling approach. The integer order modelling approach are local in nature, but, most of the fractional order modelling approach are non-local and assesses the memory effect which make itself more powerful than integer order modelling approach 10 . Due to this realistic property and their popularity to model the behavior of real systems in various disciplines, many scholars recently proposed a study of real world situations using fractional order modelling approach rather than integer order modelling approach 10,11 . Different scholars throughout the world have been formulated and analyzed mathematical modeling on various real-world situations using deterministic modelling approach or stochastic modelling approach or fractional order modelling approach see references [12][13][14] . To construct and analyze our proposed study we have reviewed some scholar studies those are relevant and constructive to our study regarding to concepts, theories, methods and methodologies. Brahim et al. 15 constructed a fractional order model to investigate the impact of harvesting on prey-predator interaction in the case of prey heard behavior. The researchers' main purpose to use fractional order model is to model the memory effect measured by the order of the fractional derivative on the mutual interactions and to investigate the impact of inner competition among predators. Teklu and Terefe 3 formulated and analyzed a susceptible, exposed, animosity-infected and treated (SEATS) compartments deterministic model on the dynamics of university students animosity towards mathematics with optimal control theory. The result concluded that applying prevention and treatment control measures simultaneously is the best strategy to minimize and possibly to eradicate the animosity-infection in the community. Mandal et al. 11 formulated and analyzed a fractional-order epidemic model with fear effect of an infectious disease with treatment control. Also fractional backward and fractional Hopf bifurcation are also analyzed and shown a role the disease control parameter, level of fear and fractional order play in the stability of equilibriums and Hopf bifurcation.
Kumar et al. 16 formulated and investigated a Caputo fractional order model on the alkali-silica reaction dynamics with six differential equations. The model analysis proved that the unique solution existence, the stability of the system using methods suitable for fractional order model approach, and illustrate the numerical simulation graphs to justify the analytical results using suitable numerical methods known as Adams-Bash forth-Moulton scheme. The study shows the significance of fractional order modelling approach in the study of chemical reactions. Erturk et al. 17 formulated a fractional order modelling approach to investigate the motion of a beam on an internally bent nanowire. Results in the model analysis shows that the fractional responses approach the classical ones as the fractional order goes to unity and additionally the fractional Euler-Lagrange equation provides a flexible model possessing more information than the classical description. Viera-Martin et al. 18 carried out a bibliographic analysis on artificial neural networks (ANNs) using fractional order modelling approach. They considered fractional order modelling approach to achieve the three objectives, such as systems stabilization, systems synchronization, and parameters training, using optimization algorithms. From the finding of the study they recommend that the Caputo fractional order modelling approach is the most utilized method for solving problems with ANNs because its initial values take the same form as the differential equations of integer-order approach. Kumar et al. 19 proposed a delay Caputo fractional order model on the oncolytic virotherapy compositing viral lytic cycle and virusspecific cytotoxic T lymphocyte (CTL) response. From the analysis of the proposed model they concluded that a fractional order modelling approach has different behaviors as compared to the integer order modelling approach and the algorithm applied in the numerical analysis part smooth and reliable for delay fractional order approach. Din et al. 20 proposed and analyzed a Caputo fractional order model on climate change. They carried out both the qualitative and numerical analysis of the proposed model and shows the result that the total spectrum lying between two integer values are achieved with more information about the complexity of the dynamics of the proposed fractional Climate Change-model. Based on the studies we have reviewed above we understand that the fractional order modelling approach could produce better solutions in the comparison of existing classical (integer order order) models, but the model analysis with fractional derivatives approach is more complicated than the classical integer order modelling approach. In this study the main and strong reason to use Caputo derivative in addition to its suitability for initial value problem, when talking about real problems, the Caputo derivative is highly useful since it allows traditional starting and boundary conditions be included in the derivation, and the derivative of a constant is zero that is not the case with the Riemann-Liouville fractional derivative 21 .
The aim of this research study is to formulate and analyze the new Caputo fractional order model with optimal control theory on the higher institutions students' anxiety towards mathematics transmission dynamics with prevention and control strategies. We are motivated with the above studies and high interest in finding the best strategies to control the higher institution students' anxiety towards mathematics using fractional order modelling approach. To the best of our knowledge, there is no scholars who studied higher institutions students anxiety towards mathematics using fractional order modelling approach so that our proposed model is unique in the given research thematic area.The remaining part of this paper is organized as follows. In section two we discussed the basic mathematical preliminaries important for the study. In section three we formulate the integer order model where each parameter and human compartment was explained. In section four we re-formulated www.nature.com/scientificreports/ and analyzed the integer order model given in section three as the Caputo fractional-order model. In section five we carried out the optimal control analysis of the Caputo fractional-order model given in section four. In section six and section seven we carried out the numerical simulations of the integer order and fractional order models respectively. In section eight we have discussed.and concluded the whole process of the research study.

Basic definitions of fractional calculus
In this section we recall some basic definitions of fractional calculus.

Definition 1
The Caputo fractional order derivative with order θ for a function h ∈ C n is defined as 22,23 Note:

Definition 2
The Caputo case fractional order integral with order α > 0 for a function h ∈ C n is defined as 22,23 Definition 3 The Mittag-Leffler function in two parameters is defined by 22,23 For θ 2 = 1 , the Mittag-Leffler function in one parameter is defined by 22

Definition 4
The constant point δ * is an equilibrium point of the Caputo-fractional model, then

Integer order model formulation
In this section, we briefly discuss the integer order model for the dynamics of higher institutions students' anxiety towards mathematics which is in the case of Caputo case fractional order derivative. In this study, to analyze the higher institutions students anxiety towards mathematics, we divide the total number of higher institution students denoted by N(t) into six distinct groups such as students who are susceptible to anxiety towards mathematics denoted by S h (t) , students who are protected from anxiety towards mathematics denoted by P h (t), students exposed to anxiety towards mathematics denoted by E h (t) , students who have anxiety towards mathematics denoted by A h (t) , students who have a permanent anxiety towards mathematics denoted by Q h (t) , and students who are recovered from anxiety towards mathematics denoted by R h (t) such that In this mathematical model formulation, since anxiety towards mathematics is not density-dependent students susceptible to anxiety towards mathematics acquire anxiety at the standard incidence rate given by Basic assumptions to develop the students anxiety towards mathematics transmission dynamics model: ε is portion of the recruited students who are entered to protected group, the susceptible group is increasing by the new recruitment portion (1 − ε) , students from protected group in which those higher institution students who are protected but who lost protection by the rate (1 − κ)ρ , κ is efficiency of protection and from recovered group who lose their temporary recovery by the rate ω , the student population is homogeneous in every group, estudents in each group are subject to natural death rate µ , student population is variable, anxiety can naturally recovered, students can recovered from anxiety after treatment measure and some students will have permanent anxiety towards mathematics. The student population is homogeneous and variable.
Using parameters in Table 1, variables in Table 2, and assumptions described above the transmission flow diagram of higher institutionse students mathematics anxiety towards mathematics is given in Fig. 1 From Fig. 1 the nonlinear ordinary differential equations governed by the model assumptions are described as follows: (2) C D θ t h(t) tends toḣ(t) as θ → 1. (3) , whereθ 1 > 0, θ 2 > 0.   The sum of all the differential equations in (Eq. 2) is dN dt = − dN.

Model derivation in Caputo fractional order derivative approach
In this section, we re-formulate the higher institutions students anxiety towards mathematics model (Eq. 9) using a Caputo fractional order derivative approach in order to observe the memory effects and gain more insights about the dynamics of higher institutions students anxiety towards mathematics. System (9) in terms of integral form is followed by substituting the value of kernel as a power-law correlation function. After applying the Caputo fractional derivative of (θ − 1), we have got Both the left and right side of Eq. (11) are the inverse operators, we have the higher institutions students anxiety towards mathematics daynamical system (Eq. 9) in the Caputo fractional order operator form The initial conditions for the fractional order derivatives (Eq. 8) are given as: Basic properties of the Cauto fractional model (Eq. 12). In the study, since we are dealing with the number of higher institutions students which cannot be negative, we have to show there is non-negative solution in the given region to the fractional order model (Eq. 12). The constructed model is mathematically analysed by proving different theorems and algebraic computation dealing with different quantitative and qualitative attributes.

Theorem 1 (Positivity of the model solutions) Each solution to the fractional order model
Proof Using similar approach stated in reference 24 , we assume by contradiction that the first equation in fractional dynamic model (Eq. 12) is not true.
Let us consider the following expression and assume it exists, Using the inverse Laplace transform on this last expression and expressed as a Mit-tag-Leffler function of one parameter, we have determine that In a similar manner one can show that for each remaing state variable and hence each solution S h (t), P h (t), E h (t), A h (t), Q h (t) , and R h (t) of the fractional order model (Eq. 12) is positive. Therefore, since we have shown the solution is bounded and positive, we can reasonably say that the model is both biologically and mathematically meaningful.
Proof By adding all the fractional order derivatives in Eq. (12) we have determined as . Simplifying the result to determine the following inequality L( And solving using the definition of inverse Laplace transformation and using definition (Eq. 5) we have determined the expression given by for each t ≥ 0. Therefore, the total number of students population N(t) is bounded and hence the proposed fractional order model (Eq. 12) is both mathematically and biologically meaningful.
Qualitative analysis of the fractional order model (Eq. 12). Anxiety-free equilibrium point of the model (Eq. 12). For the higher institutions students anxiety-free equilibrium point of the fractional order model (Eq. 12), we have calculated from CD θ Effective reproduction number of the model (Eq. 12). In this manuscript the next generation matrix operator approach by Van den Driesch and Warmouth as in 25 is used to derive the effective reproduction number R θ eff for the dynamics of the higher institutions students anxiety towards mathematics. The higher institutions students anxiety towards mathematics model (Eq. 12) effective reproduction number denoted by R θ eff measures the average number of the higher institutions anxiety infected students generated by one anxiety infected student in a considered community of students when some controlling strategies are in place, like anxiety protection, or/and anxiety treatment. The students anxiety infection effective reproduction number denoted by the symbol R θ eff is the dominant (largest) eigenvalue (spectral radius) of the matrix where F i is the rate at which newly anxiety infected students appear in compartment i , ν i is the transfer of students who have anxiety towards mathematics from existing compartment i to another compartment and E 0 A is the students anxiety-free equilibrium point. Let and after a detailed computations, we have got the transmission matrix as. and the transition matrix as The next generation matrix is of the form , P 0 = ε θ � θ µ θ +(1−κ θ )ρ θ and D 1 = µ θ + ψ θ µ θ + δ θ + γ θ . Then the spectral radius of the matrix FV −1 which is called the higher institutions anxiety infected students model (Eq. 12) effective reproduction number given by Local stability of anxiety-free equilibrium point.
Theorem 3 Given a fractional-order system of differential equation D θ 0 y(t) = f y , 0 < θ ≤ 1. Let y 0 be an equilibrium point of the dynamical system, and let B = D(f y 0 ) be the Jacobian matrix of f evaluated at y 0 . Then y 0 is locally asymptotically stable if and only if arg( i ) > θπ 2 , for each eigenvalue of the matrix B 22 .
Note: The proof of Theorem 2 illustrates that in case 0 < θ < 1 , the stability region of the fractional-order model (Eq. 12) increases as compared to the integer order model.

Theorem 4 The anxiety-free equilibrium point E 0
A of the fractional order model (Eq. 12) is locally asymptotically Proof The local stability of the anxiety-free equilibrium point of the model (Eq. 12) at the anxiety-free equilib- , 0, 0, 0 has been studied using the criteria stated in Theorem 3.
The Jacobian matrix of the dynamical system given in (Eq. 12) at the anxiety-free equilibrium point is given by The computations for eigenvalues of the Jacobian matrix J E 0 A gives and where We can justify that a 0 , a 1 ,a 2 and a 3 are all positive implies all the three eigenvalues of Eq. (14) have negative real parts if R θ eff < 1. Then since all the eigenvalues of the Jacobian matrix J E 0 A have negative values using the criteria stated in Theorem 3 above arg( i ) > θπ 2 for 0 < θ ≤ 1 and each i = 1, 2, 3, 4, 5, and 6. Thus, the anxietyfree equilibrium point of the model (Eq. 12) is locally asymptotically stable whenever R θ eff < 1.
Existence of anxiety endemic equilibrium point (s) and bifurcation analysis. The anxiety endemic equilibrium point occurs when the anxiety exist in the higher institution student population. Anxiety endemic equilibrium point of the Caputo fractional order derivative model (Eq. 12) is obtained by making right hand side of all equations of system (Eq. 12) equal to zero.
be the students anxiety standared inci- www.nature.com/scientificreports/ dence rate ("force of anxiety infection") at the anxiety endemic equilibrium point. To find equilibrium point(s) for which anxiety infection is endemic in the population, the model equations given in Eq. (12) are solved in  and Since all the model parameters are non-negative one can be seen from Eqs. (17) and (19) that a 2 > 0 . Moreover, a 0 > 0 whenever R θ eff < 1 . Thus, the number of possible positive real roots the degree two polynomial (Eq. 18) can have depends on the sign of a 1 . We analyzed it using the Descartes' rule of signs on the quadratic polynomial f (x) = a 2 x 2 + a 1 x + a 0 (with x = * n ).

Theorem 5
The higher institutions students anxiety towards mathematics fractional order model (Eq. 12).
(a). has a unique endemic equilibrium if R θ eff > 1 and either of the following holds, (i).a 1 > 0, (ii). a 1 < 0, (b). Could have two endemic equilibrium if R θ eff < 1 whenever a 1 has negative sign . Theorem 5 condition (b) suggests the possibility of existence of multiple endemic equilibriums whenever R θ eff < 1 (which is typically shows the existence of the phenomenon of backward bifurcation see references 23,[26][27][28][29][30][31][32][33] . The phenomenon of backward bifurcation is characterized by the co-existence of a stable anexiety-free equilibrium and a stable anxiety endemic equilibrium whenever the effective reproduction number of the model is less than unity. .

Theorem 6 Suppose E *
A be the anxiety endemic equilibrium point of the model (Eq. 12) given by Eq. (16) and stated in Theorem 5 then it is locally asymptotically stable. Proof To prove the local asymptotic stability of the anxiety endemic equilibrium poin we can apply the fractional odrder Routh Hurwitz stability criteria stated in reference 34 , it is enough to show that all eigenvalues of the following Jacobian matrix satisfy the Matignon condition 35 stated in Theorem 3.
Global asymptotic stability of the anxiety endemic equilibrium point.

Theorem 7 Let 0 < θ < 1 be the order of the fractional order model given in (Eq. 12) then the unique anxiety endemic equilibrium point E *
A whenever R θ eff > 1 is globally asymptotically stable.
Proof Based on the Lyapnove function method stated in 23,36 , we can construct the following function Since H(t) is contionous the Caputo fractional order derivative of H(t) is computed as.
Then after solving this expression we have determined the following result and a 5 = µ θ + δ θ + γ θ µ θ + ψ θ µ µ θ + ω θ * Hence, by 22,23 the function H is a Lyapnov function on the feasible domain and the largest set in this feasible domain satisfies the condition Therefore, the anxiety endemic equilibrium point is globally asymptotically stable if all these conditions and R θ eff > 1 holds.

Optimal control problem
In this section, we extend the fractional order model (Eq. 12) by introducing two time-dependent controlling strategies, where c 1 (t) represents efforts to prevent students from anxiety infection and help reduce anxiety contact rates and c 2 (t) represents the intensity of educational treatment of mathematics anxious students to increase recovery from anxiety where 0 ≤ c 1 (t), c 2 (t) ≤ 1 . It is assumed that the exposed population in susceptible students is reduced by the factor (1− c 1 (t) ) due to the protection measures taken. Similarly, the anxiety infected students is reduced by the factor (1− c 2 (t) ) due to the educational treatment by experts. Hence, the control theory dynamical system of the Eq. (12) becomes: To minimize the number of mathematics anxious students in the community we construct the objective function given by The control problem involves a situation in which the number of mathematics anxiety infected students and the cost of applying preventions and treatments controls u 1 (t) and u(t) are minimized subject to the system (29). Where T is the final time, the coefficients ℧ 1 and ℧ 2 are positive weight constants and B 1 2 and B 2 2 are the measure of relative costs of interventions associated with the controls c 1 and c 2 , respectively, and also balances the units of integrand. The objective is to find the optimal values C * = c * 1 , c * 2 of the controls C = (c 1 , c 2 ) such that the associated state trajectories h are solution of the system (Eq. 22) in the intervention time interval [0, T] with initial the given conditions and minimize the objective functional. In the cost functional, the term ℧ 1 E h refer to the cost related to anxiety exposed students and the term ℧ 2 A h refer to the cost related to anxiety infected class.

Theorem 8 (Existence of optimal solution)
There exists an optimal control C * = c * 1 , c * 2 in C and a corresponding solution vector Note: We utilize Pontryagin's maximal principle stated in 3,37 to determine the prerequisites for the optimal control model (Eq. 22). For the optimal control problem (Eq. 22) we define the Hamiltonian (H) function by www.nature.com/scientificreports/ where 1 (t), 2 (t), 3 (t), 4 (t), 5 (t) and 6 (t) are the co-state variables or adjoint variables. Using the same method stated in 37 for fractional order model approach we have determined the following: The transversality conditions are On the interior of the control set, where 0 < c i < 1 for i = 1, 2 we have the following equations Then solving for for c 1 and c 2 gives us
The optimality conditions is obtained by differentiating Hamiltonian H with respect to the control variables c 1 and c 2 : ∂H ∂c 1 = 0 =, ∂H ∂c 2 = 0.

Numerical simulations for the deterministic model (Eq. 9)
To illustrate the numerical results of the integer order model (Eq. 9), we consider the fixed parameters values some of them are stated in reference 3 and some of them are assumed and given by = 100, µ = 0.5, ψ = 0.04,ε = 0.4, κ = 0.8, γ = 0.01, ω = 0.03, ρ = 0.2,δ = 0.3, σ = 0.04, ϕ = 1.3 , and β = 1.4 where some of these values are taken from the animosity study in reference 3 . In this section we execute a numerical simulation of the higher institution students anxiety towards mathematics model (Eq. 3) to justify the analytical results we performed in "Model derivation in Caputo fractional order derivative approach" using Matlab standard ordinary diferential equations (ODEs) solver function ode 45 .  Fig. 3 investigated the effect of anxiety transmission rate β on the anxiety effective reproduction number R eff . The figure expresses that when the value of the anxiety transmission rate β increases, the anxiety effective reproduction number increases, and whenever the value of β < 0.601 implies that R eff < 1. Therefore, the responsible body shall concentrate on minimizing the value of the anxiety transmission rate β to prevent and control anxiety transmission in the higher institution community.  Fig. 4 illustrated that the effect of anxiety towards mathematics treatment rate γ on the anxiety effective reproduction number R eff . The figure shows that when the value of the anxiety treatment rate γ increases, the anxiety effective reproduction number deacreases, and whenever the value of γ > 0.832 implies that R eff < 1. Therefore, the concerned body shall concentrate on maximizing the values of anxiety treatment rate γ to minimize anxiety spreading throughout the student community.

Numerical simulations of the fractional order model (Eq. 12)
In this section the numerical simulation of the Caputo fractional order model (Eq. 12) is simulated by using fractional Euler's forward method. To illustrate the numerical results of the fractional order model (Eq. 12), we considerd two possibilities for the order of the derivative as θ = 0.5 and θ = 1 and we have taken the fixed parameters used in the simulation of the integer order model (Eq. 9) with values = 100, µ = 0.5, ψ = 0.04,ε = 0.4, κ = 0.8, γ = 0.01, ω = 0.03, ρ = 0.2,δ = 0.3, σ = 0.04, ϕ = 1.3 , and β = 1.4 where some of these values are related to animosity study in 3 . In this section we execute a numerical simulation of the higher institution students anxiety towards mathematics model (Eq. 12) to justify the analytical results we performed in "Model derivation in Caputo fractional order derivative approach" using Euler forward method and writing a Matlab code for the fractional order diferential equations (FODEs) given in (Eq. 12). In order to observe the effects that the parameter θ has on the dynamics of the fractional-order model (Eq. 12), we include several numerical simulations varying the value of this parameter.

Effect of memory on higher institutions students anxiety infection.
Here, we simulate the order of derivative (memory) effects ( θ ) on the number of anxiety exposed, infected and recovered higher institutions students. In the simulations given in Figs. 7, 8, 9 and 10, we compare the number of anxiety exposed, infected and recovered higher institutions students for memory-less(integer order) model that is when ( θ = 1 ) and model with memory that is when ( θ = 0.5 ). We illustrated the effects of memory(order of derivatives)(θ ) on the number of anxiety exposed, infected and recovered higher institutions students from Figs. 7, 8, 9 and 10. As we observe in Fig. 7, the number of anxiety infectious students (A h ) is larger in the system with out memory ( θ = 1 ) as compared with system with memory ( θ = 0.5 ). Similarly from Figs. 8, 9 and 10 we observed that the number of anxiety exposed (E h ), permanently anxitious (Q h ) and anexiety recovered (R h ) are larger in system with out memory ( θ = 1 ) compared with systems with memory ( θ = 0.5 ). In this section computation of the model (Eq. 12) effective reproduction number is determined as R θ eff = 2.83 which shows the persistence of mathematics anxiety in the student community.
Effect of optimal control strategies on higher institutions students mathematics anxiety. In this sub-section we have carried out the numerical simulation of the state variables in the optimal control problem given in Eq. (21) to investigate the effect of controlling strategies using Euler forward method where the order of the derivative is assumed to be θ = 0.5 and some of the fixed parameters values given in the simulation of the deterministic model (Eq. 9) as = 100, µ = 0.5, ψ = 0.04,ε = 0.4, κ = 0.8, γ = 0.01, ω = 0.03, ρ = 0.2,δ = 0.3, σ = 0.04, ϕ = 1.3 , and β = 1.4 where some of these values are related to the animosity study in 3 . Numerical simulations illustrated by the following figures show the significance of the control strategies to takle the transmission dynamics of mathematics anxiety in the students community.
Simulation illustreated by Fig. 11 shows that when the protection strategy increases then the number of students who are susceptible to mathematics anxiety (S h ) decreases while the number of students who are protected against mathematics anxiety (P h ) increases.
Simulation illustreated by Fig. 12 shows that when the rates of protection and treatment strategies increases then the number of anxiety exposed students (E h ) decreases.
Simulation illustreated by Fig. 13 shows that when the rates of protection and treatment strategies increases then the number of anxiety infected students (A h ) decreases. Based on the numerical simulation results we ovserbed that applying both protection and treatment controlling strategies simultaneously has a fundamental effect on the transmission dynamics of mathematics anxiety throughout the higher institutions students community.

Discussion and conclusion
In this study, we formulated and analyzed a deterministic mathematical model and reformulated it as a fractional order mathematical model on the higher institutions students' anxiety towards mathematics with prevention and treatment mechanisms. In "Introduction" we have introduced the basic concepts students anxiety towards mathematics and some basic background for this study. In "Basic definitions of fractional calculus" we have recalled some basic definitions of fractional calculus which are fundamental for this resaerch study. In "Integer order model formulation" we have formulated and briefly discuss the integer order model of the transmission dynamics of higher institution students anxiety towards mathematics using a system of ordinary differential equations by dividing the total number of higher institution students population into six distinct groups. In "Model derivation in Caputo fractional order derivative approach", we have re-formulated the integer order model given in "Integer order model formulation" into a Caputo fractional order approach and analyzed the qualitative behaviors of the model such as; non-negativity of future solutions of the model, boundedness of the dynamical system, existence of anxiety-free equilibrium point, existence of anxiety effective reproduction number using next generation matrix approach, existence of anxiety endemic equilibriums, local stability analyses of anxiety-free and anxiety endemic equilibrium points using Routh-Hurwiz and Matignon's stability criteria for fractional oreder model. Backward bifurcation in fractional order approach is estabilished. In "Optimal control problem" an optimal control problem for the fractional order dynamical system counterpart is re-formulated www.nature.com/scientificreports/ and investigated. In "Numerical simulations for the deterministic model (Eq. 9)" numerical solutions for the deterministic model is carrieds out. In "Numerical simulations of the fractional order model (Eq. 12)" numerical simulation for the fractional order approach model including its counterpart optimal control problem has been performed and used to verify the qualitative (theoretical) analyses of the model. The reason for we considered a fractional order approach instead of its integer order counterpart is that the fractional order differential equation approach is a generalization of integer order differential equation. One can argue that a fractional order approach is more suitable approach than integer order approach for modelling any complex adaptive systems in different diciplines of study. We have carried out numerical simulations for both the deterministic and fractional order approaches using MATLAB programming codes with ODE45 (the fourth order Runge-Kutta) approach to the simulation of deterministic model and with Euler forward finite difference approach for the fractional order model. From the numerical simulation results we illustrated that some parameters changes have high impacts on the anxiety effective reproduction number R eff of the model and determined the result R eff = 5.41 which shows the persistence of high anxiety in the student community, the behavior of the deterministic model solusions and effects of some influential parameters like transmission rate, protection rate and treatment rate on the model solutions, we computed the effective reproduction number of the model (Eq. 12) as R θ eff = 2.83, investigated the impact of memory on the number of anexiety exposed, anxiety infectious, permanently anxious and recovered form anxiety students. Also from the numerical simulation results of the optimal control problem we ovserbed that applying both the protection and treatment strategies is the most effective approach to tackle the mathematics anxiety in the students community. In general, our fractional order model numerical simulation result shows that memory has great influence on anxiety infection transmission. Even though, we suggested that the fractional  Time in years Anxiety exposed students E h Effect of control strategies on anxiety exposed tudents c1=0, c2=0 c1=0.5,c2=0.6 Figure 12. Effect of protection and treatment strategies on the number of anxiety exposed students ( E h ). www.nature.com/scientificreports/ order modelling approach could produce better solutions in the comparison of existing classical models, we strongly believe that this study analyses can further modified by potential researchers. However, we understand that the model analysis with fractional derivatives approach is more complicated than the classical deterministic modelling approach. To the best of my knowledge this is the first paper on higher institutions students anxiety towards mathematics in fractional order modelling approach. Eventually, since protection rate and tereatment rate have fundamental impacts to minimize the transmission dynamics of higher institution students anxiety towards mathematics we recommend for stakeholders to concentrate on maximization of both protection and treatment measures to tackle higher institutions students' anxiety towards mathematics. Finally, since this research study is not exhaustive any potential researcher can modified this research study by incorporating additional concepts such as the stochastic approach, age structure of students, effects of teaching aid materials, roles of parents, and fitting the model with real situations data.

Data availability
Data used to support the findings of this study are included in the article.